Inflationary theory of branching morphogenesis in the mouse salivary gland

The mechanisms that regulate the patterning of branched epithelia remain a subject of long-standing debate. Recently, it has been proposed that the statistical organization of multiple ductal tissues can be explained through a local self-organizing principle based on the branching-annihilating random walk (BARW) in which proliferating tips drive a process of ductal elongation and stochastic bifurcation that terminates when tips encounter maturing ducts. Here, applied to mouse salivary gland, we show the BARW model struggles to explain the large-scale organization of tissue. Instead, we propose that the gland develops as a tip-driven branching-delayed random walk (BDRW). In this framework, a generalization of the BARW, tips inhibited through steric interaction with proximate ducts may continue their branching program as constraints become alleviated through the persistent expansion of the surrounding tissue. This inflationary BDRW model presents a general paradigm for branching morphogenesis when the ductal epithelium grows cooperatively with the domain into which it expands.


Supplementary Information
Inflationary theory of branching morphogenesis in the mouse salivary gland Here we develop in more detail the computational and theoretical framework used to study the branching dynamics of the salivary gland ductal epithelium. First, we define the agent-based implementation of the stochastic IBDRW model, and how the results of stochastic simulations were compared and fit to the experimental data. We then turn to the construction of the mean-field approximation of the IBDRW model, which we use to gain analytic insight into the growth dynamics of the branching process. Finally, we discuss in brief the potential application of the IBDRW to other branching contexts.

S1. STOCHASTIC SIMULATIONS OF THE IBDRW MODEL USING AN AGENT-BASED REPRESENTATION
To investigate the capacity of the BARW and its generalization, the IBDRW model, to predict the statistical distributions of the salivary gland network, we implemented a computational agent-based model of the stochastic branching process. The IBDRW comprises a set of processes illustrated in Fig. 4a, which include stochastic tip propagation, duplication (tip bifurcation), delay and reactivation, which all take place on the background of a global, constant, inflation of the domain on which the dynamics takes place. From the biological perspective, this means that, during development, ducts must expand through cell proliferation to match the tissue expansion. Therefore, in the context of the salivary gland, the distance between any two distinct points on the network necessarily increases exponentially over time (see Figs. 1 and 2f).
To simulate the full stochastic dynamics of the IBDRW model, we considered an off-lattice, discrete time, agentbased implementation. The dynamics of the conventional BARW model follows the same program but in the absence of domain expansion. To capture the dynamics, three particle species were considered: active tips (endbuds), immobile tracer particles (representing the locus of trailing ducts), and delayed tips (endbuds partially inhibited by local steric constraints). Active tips were subject to two concurrent Poisson processes: tip branching and propagation.
Formally, at each time-step of fixed duration ∆t, we determined the number, a(t), of active particles at time t. We then constructed a list of uniformly distributed random numbers r i on the interval [0, 1], with i = 1, ..., a(t), and cycled through all active particles in random order, allowing them to either branch or "hop" in the manner of a persistent random walk. A particle j underwent branching if the branch probability r j < p b = (1 − exp(−R branch ∆t)) (with R branch the constant branch rate), or otherwise hopped. Branching was implemented non-locally such that, when a particle branched, it was removed and replaced by two active tips at new locations.
Following the propagation of an active tip, an inactive particle was left behind, representing the trail of elongating ducts. Active particles became delayed when they came within an inhibition radius, ρ a , of a tracer particle. We also applied a "sensing cone" such that an active particle was only inhibited by tracers laying within a spherical sector around its direction of motion. This had two purposes: first, this prevented the (non-physical) immediate inhibition of active tips due to their interaction with their own proximate trailing duct; and, second, this ensured that inhibition did not immediately occur between two newly-born tips, but was rather due to direct tip-duct interactions. (More generally, in systems where inhibition occurs in response to biochemical signals secreted from maturing ducts, tips may sense chemical diffusion gradients coming from any direction. However, here, we focused on a model in which tips are inhibited by steric influences alone.) The ductal network of the salivary gland develops as a three-dimensional tissue. We therefore considered the tip growth dynamics as a 3d persistent random walk [1], with zero translational diffusion but non-zero rotational diffusion constant D r . Given the 3d position r i (t) = (x i , y i , z i ) of an active tip at time t, with a director oriented according to the polar and azimuthal angles (θ, ϕ) in spherical coordinates, an elongation event implied a translation of the active tip from position r i (t) to position r i (t + ∆t) = r i (t) + ∆L(sin θ ′ cos ϕ ′ , sin θ ′ sin ϕ ′ , cos θ ′ ), where the step length, ∆L = V 0 ∆t, corresponds to the displacement of a tip, moving at a speed V 0 , during the time-step ∆t, and (θ ′ , ϕ ′ ) = (θ, ϕ) + 2D r ∆t arccos(1 − 2u 1 ) π − 1/2, u 2 .
Here, the second term on the right-hand side corresponds to a uniformly distributed noise on a 3d sphere, with u 1 ∈ U ([0, 1[) and u 2 ∈ U ([−1, 1]), with amplitude √ 2D r ∆t. The angular diffusion constant D r encodes the persistence time of ductal growth. In the absence of noise, D r = 0, ducts grow along a straight line, exhibiting infinite persistence. However, for D r > 0, ducts diffuse with τ = D −1 r denoting the persistence time and ℓ = V 0 τ the persistence length. The former corresponds to the typical time it takes for the direction of motion of a tip to become uncorrelated, providing a natural time scale for the growth of the spatial branching pattern. Varying the time scale τ and branching rate R branch accordingly results in networks with identical topology.
During a branching event, two new tips were generated. The parent duct and the two offspring branches were considered to be co-planar, while the angle of rotation or roll of this plane around the axis formed by the parent duct was estimated empirically (see Supplementary Figs. 6a-6e), as was the angle formed between the two offspring branches (see below and Supplementary Fig. 6d).
Exponential domain growth was implemented as follows: at each time step, the positions r i (t) of every particle (tip or duct) in the system, were scaled by a factor (1 + R exp ∆t), which corresponds to a long-term exponential growth of the original domain size in all three dimensions by a factor ω(t) = exp(R exp t), with expansion rate R exp . To account for reactivation of tips, we kept track of all delayed tips in the system and checked, on every iteration and after the expansion of the domain, whether they were still within the inhibition radius of a duct "particle". If this was not the case, the tip was reactivated, allowing it to undergo branching and elongation events. For every realization of the process, we recorded the positions of all particles, their identity and connectivity in the form of an edge list. This allowed the reconstruction of the branching tree to study their spatial and topological properties and compare them with the experimental observations.

A. Initialization, boundaries and tip termination criteria
Depending on the context, stochastic simulations were initialized either by placing a single duct at the origin of the expanding domain or by seeding the system with one of the measured ductal tree configurations at E14.5. In the case of the former, the initial duct had only one active tip, with the other end considered as the root of the branching tree. In the latter, statistical distributions were obtained from the average of many stochastic simulations with 1/3 seeded with the measured E14.5 configuration of sample 1, 1/3 from sample 2 and 1/3 from sample 3 (see Supplementary Fig. 3). In all cases, at the initial time-point, all tips were taken as active except for the root of the tree. Stochastic simulations were carried out in an unconfined (open) domain, which allowed tips to propagate outwards without boundary constraint until the maximum simulation time t max was reached. This choice was consistent with the hypothesis that, during embryonic development, the salivary gland expands and branches freely until there is a global transition to a phase of acini formation, commencing after E16.5, that terminates fully the branching process [2].
B. Estimation of model parameters from the experimental data As discussed earlier, the persistence time (τ = 1/D r ) provides a time scale for the growth of the branching pattern. With this in mind, the IBDRW has only three key parameters: the ratio of branching rate to rotational diffusion r branch = R branch /D r , the inhibition radius ρ a , and the ratio of the gland expansion rate to rotational diffusion r exp = R exp /D r . Before estimating the rates, we first extracted the characteristic length and time scale of the tip growth process by estimating the level of persistence of ductal elongation. For this, we measured the duct tortuosity [3] at E14.5. Given the duct length L d between two branching points, and the Euclidean distance ℓ d between the two, we defined the duct tortuosity as the relative deviation ∆ = (L d − ℓ d )/L d ≈ 0.11. This implied that, on average, ducts were 11% longer than the Euclidean distance between consecutive branching points. For small deviations, the value of the rotational diffusion constant D r was approximated by ∆, so that the factor √ 2D r ≈ π/7 for the salivary gland, with a characteristic time scale of persistence τ = 9.1 h.
To estimate r branch , we first measured the average level reached by each endbud in E14.5 and E18.5 trees. Sampled across n = 3 replicates at E14.5, endbuds were, on average, positioned at level 9 ± 1 (average±s.d.), while at E18.5 this increased to level 14 ± 1 (average±s.d.), resulting in a typical difference of about 5 ± 2 (average±s.d.) generations of branching between the two time-points. This indicated that a typical active endbud at E14.5 undergoes around 5 subsequent branching events in the course of 4 days. This yielded a tip branching rate of around R branch = 0.05 ± 0.02 h −1 (average±s.d.) or 1 tip bifurcation event every 19 hours. Rescaling by τ −1 , we therefore have r branch = 0.45.
We then proceeded to estimate the duct elongation speed, V 0 . For this, we considered the average duct length between consecutive bifurcation events at generations 15 or above. This choice was made to minimize the effect of tissue and duct expansion, which would otherwise cause an effective reduction of the estimate. This resulted in an average duct length of 20 ± 2µm (average±s.d.), or taking into account the time between branching events, an elongation speed of V 0 = 1µm/h. An estimate of the expansion rate R exp was made by measuring the typical diameter of the SMG at different time points (Fig. 2). By fitting exponential growth curves, proportional to exp(R exp t), we found that the growth rate R exp for the whole SMG was given by R exp ≈ 0.014 h −1 , 95% CI [0.011, 0.017] (n = 3 repeats for each time point). When focused on the first lobe, we found a reduced rate of R exp ≈ 0.0135 h −1 , 95% CI [0.0066, 0.020]. From this, we considered the expansion for the SMG to be R exp ≈ 0.01 h −1 , or rescaling by τ −1 , r exp ≈ 0.09.
To determine the inhibition radius, a rough estimate was made based on the distance between endbuds and the nearest duct, in a conic section in their direction of propagation. This provided an estimated value of ρ a ≈ 60 ± 34 microns (average±s.d.). An alternative estimate, based on the nearest-neighbour distance between acini at the E18.5 time point, resulted in the value ρ a ≈ 36 ± 26 µm (average±s.d.). Note that, in both cases, the range of values of the inhibition radius are broad. This is because its estimate is subject to sources of technical variability. In particular, measurements of tip-to-duct and tip-to-tip distances were carried out by measuring the distance from the center of a duct to the center of an acinus, as this is the quantity that enters in the agent-based simulations. However, the actual size of a duct and acinus may vary from duct to duct, which will affect the estimate of the range of inhibition. As such, the inhibition radius of the model may not translate precisely to the measured estimate. Moreover, within the IBDRW framework, we assumed that tip delay responds only to tip-duct interaction, i.e., volume exclusion. However, in the context of the real biological system, we note that there might be biochemical cues and more complex interactions that contribute to an effective inhibition radius. With these caveats in mind, based on the estimates above, we considered an rough estimate of the inhibition radius of around 50 microns, a value that was then used as a guide to define the reasonable fitting range in the stochastic simulations.
The total developmental time (which translates to the maximum simulation time) was set to 84 hours, or 3.5 days. This assumes that differentiation of tip cells commences at around E18, so that branching is reduced and by E18.5 all endbuds have differentiated into acini. Larger values of the developmental time resulted in larger and more persistent subtrees (see Supplementary Fig. 9e).
The angles formed by offspring branches during a tip bifurcation were estimated from the empirical data. The (periodic) roll angle φ of the plane formed by the two offspring branches followed a peaked distribution that we approximated by a Laplace distribution with mean µ = 90 • and scale parameter b = 21 • (Figs. 6a-6e). Empirically, the angle formed by two offspring branches at E15.5 was well-fit by a truncated normal distribution, with mean 90 • and standard deviation 20 • , which increased slightly at E18.5 ( Supplementary Fig. 6d). Small changes to these parameters had no significant effect on the quantitative properties of the networks. However, as observed in other contexts [3,4], large angular deviations do affect the organization, mainly because small branching angles may lead to immediate termination of one or both offspring tips. The quantitative properties remained, to a large extent, unaffected as long as the branching angle and sensing cone angles did not cause immediate (and non-physical) termination events after branching. To find a best fit of the IBDRW to the experimental data, we generated stochastic simulations of the model for the parameter values defined above, taking as an initial condition one of the three measured rudimentary trees at E14.5. The inhibition radius was varied in the range ρ a ∈ [30, 70]µm and the expansion rate in the range r exp ∈ [0.003, 0.014]. For each set of parameters, we computed the endbud fraction q i and distributions of subtree size and persistence from N = 10 3 realizations. By comparing the results of stochastic simulations with measurements from the experimental data using an R 2 goodness-of-fit estimate (see the following section), we found the parameter values of r exp and ρ a that provided the best fit to the data. The summary of parameters are shown in Table S1. Alongside fits to the full branching statistics, further predictions were made of the variation in the average duct length with branch level, the duct length distribution and density fluctuations, and the results compared with experiment.
To test whether salivary gland development was consistent with a conventional BARW model, we considered the same parameter estimations defined above for r branch and ρ a , supplemented by the branching angles and developmental time, and considering the ductal networks at E14.5 as the initial conditions of the numerical simulations (see Table S1). However, in this case, we set r exp = 0. As noted in the main text, this model failed to recapitulate the spatial and topological properties of the true ductal networks ( Supplementary Fig. 7). In particular, this model predicted much smaller subtrees, and shorter duct lengths compared to the E18.5 data.
Finally, as emphasized in the main text, the inability of the conventional BARW to predict the growth dynamics of the salivary gland was not surprising. The conventional BARW model predicts two essential features: a linear-like expansion of the network and the self-organization of tip growth to the boundary. By contrast, the salivary gland shows an exponential-like growth characteristic, with proliferatively active tips observed throughout the network during development.

C. Goodness of fit
Fits of the stochastic simulations of the IBDRW model to the experimental data were performed using the coefficient of determination R 2 . Given a set of n points to fit in the data (x 1 , x 2 , ..., x n ), whose average value is given byx, and a corresponding set of theoretical predictions (y 1 , y 2 , ..., y n ), and defining R 2 takes the value of unity if the prediction fits the data perfectly. In light of the exponential-like decay of many of our observations, such as the duct length distributions (Fig. 2a) and subtree statistics, we computed R 2 log , defined as the coefficient of determination with the log-transformed inputs x i and y i . This gives more weight to those values in the tails of the data, thus providing a better measure of the goodness-of-fit.

D. Sensitivity analysis
To test how strongly the model predictions are influenced by changes to each of the input parameters, we ran stochastic simulations of the full stochastic IBDRW model, choosing a range of parameter values above and below those reference values shown in Table S1, based on direct measurements or fits. The results of this sensitivity analysis (Supplementary Fig. 9) showed that all tested parameters, with the exception of the duct rotational diffusion D r , had an important effect on the subtree sizes distribution. Variations of the branching rate R branch had a noticeable effect on all statistics ( Supplementary Fig. 9a), while variations of the expansion rate had the most effect on the duct length distribution ( Supplementary Fig. 9b). Variations in the simulation time ( Supplementary Fig. 9d) led to changes in the overall tree sizes and persistence. Similarly to the branching rate, variations of the annihilation radius also affected all statistics. On the other hand, changes to the rotational diffusion had little effect on the spatial and topological measures ( Supplementary Fig. 9e).

E. Zero-dimensional stochastic branching and termination model
As an alternative hypothesis to explore the branching dynamics, we considered a zero-dimensional model in which the chance of branching and termination was assumed to be autonomous, regulated at the level of individual tips, and independent of spatial location within the tissue and local duct density [5,6]. In this scenario, the statistical organization of the salivary gland branching tree is characterised fully by the probability q i that a tip terminates at a given level, here referred to as endbud fraction. Mathematically, this model corresponds to a stochastic birthdeath (BD) type process, where the branching probability depends on the generation in the tree. To test this model, we implemented numerically a discrete model in which particles a, representing active tips, were allowed to branch a −→ a + a with probability 1 − q i or become inactive (s), a −→ s with probability q i , where q i depends on the level i and was extracted empirically from the experimental data at E18.5, as depicted in Supplementary Fig. 7a.
Numerically, stochastic simulations were performed by initialising the system with a single particle at level 1, and on every iteration an active tip, chosen at random, was allowed to branch or terminate with a probability that depended only on its branch level. Simulations terminated once all tips in the tree became inactive. We compared the empirical observations with the predictions of the model from N = 10 3 realizations (see BD in Figs. 7a and 7d). These results show that a simple birth-death model fails to recapitulate the statistical properties of the salivary gland branching tree, predicting drastically smaller but more persistent trees than found experimentally. These observations suggest that spatial, extrinsic, constraints are key regulators of the topological organization of the salivary gland epithelium.

S2. MEAN-FIELD THEORY
The agent-based implementation of the IBDRW model allowed us to quantify the statistical properties and dynamics of the model and compare the results with experimental measurements from the developing salivary gland network. However, such an implementation does not provide analytical insights into the dynamics of the growth process, nor can it identify which parameter combinations are key to the dynamics. To gain further insight into the dynamics of the IBDRW, we turned to an effective hydrodynamic or mean-field theory of the stochastic branching process. The mean-field theory can be motivated as a coarse-grained description of the dynamics or, more formally, obtained as the minimum action of the non-equilibrium field theory of the stochastic IBDRW process [7]. Although, the mean-field approximation is not capable of predicting the statistical organization of the network, it can define how the local average density of active and inactive endbuds, as well as ducts, evolve as a function of developmental time.
A. Mean-field theory of the BARW model To develop the mean-field theory of the IBDRW model, we start by first defining the set of "reactions" for the composite components of the network in the conventional BARW. Formally, active tips (a) are modeled as point-like particles that undergo a persistent random walk. As an active particle moves, it leaves behind a trail of inactive ducts (i), that constitute the ductal network. When an active tip becomes within an inhibition radius of an active tip or duct, it becomes inactive. This process can be represented as the set of local reactions [ where R branch , R duct , and R delay correspond to the branching, duct creation, and inhibition rates, respectively. In this framework, the inhibition radius need not be incorporated explicitly. Instead, we consider that all reactions occur locally, and that the effect of the annihilation radius is subsumed into the effective inhibition rate, R delay . Here, we consider only dichotomous branching in which an active tip branches into precisely two new active tips. This could be generalised by considering a distribution in the number of offspring, P (n), with n ∈ N, so that a −→ na. However, it has been shown that the long-term dynamics of such systems is always dominated by dichotomous branching [8].
Approximating the motion of active particles by a diffusive process with diffusion constant D, the set of reactions (S5) may be approximated, at the mean-field level, by a set of partial differential equations (PDEs) describing the spatio-temporal evolution of the space-and time-dependent densities, A(x, t) = [a] and I(x, t) = [i], of active and inactive (i.e., duct) particles [3], respectively: In the framework of the BARW model (S6), tip-driven growth is constrained by steric effects on a fixed domain, so that temporal changes in the ductal network occur only in response to the activity of tips. Active tips submerged within the trailing ductal network become rapidly terminated so that growth becomes driven by proliferation of active tips that self-organize as a soliton pulse at the expanding front of the epithelial network. As this pulse propagates, it leaves in its wake a spatially uniform distribution of inactive ductal tissue. As noted in the main text, more complex dynamics, such as a persistent random walk (or active Brownian motion), may be introduced at the mean-field level [9]. However, as these refinements enter as gradient terms, they would not effect the long-term steady-state dynamics of the expanding network.
B. Mean-field theory of the IBDRW In the IBDRW model, branching of the ductal epithelium takes place on the background of an expanding domain ( Figs. 1 and 2f). The expansion of the domain has three effects on network development: 1) the promotion of a constant linear-expansion of ducts to match precisely the expansion of the tissue and preserve the connectivity of the network; and 2) a local decrease of epithelial density, which leads to 3) the reactivation of delayed tips due the increase in available space. The second effect is accommodated at the mean-field level when considering the dynamics on an expanding domain. Persistent duct expansion and tip reactivation may be incorporated directly into the reaction-diffusion scheme by introducing a third particle species s corresponding to delayed tips that are generated upon inhibition of active tips and reactivated at rate constant R react . In this case, the set of reactions for the IBDRW model takes the form The fact that ducts co-expand with the domain means that the duct expansion rate R elong coincides precisely with the linear expansion rate of the system, i.e., R elong = R exp . Moreover, since the reactivation of endbuds follows directly from the expansion of tissue as new space is opened, we may set R react = R exp .
Before incorporating the background expansion, at the mean-field level, the set of reactions (S7) may be cast as a set of kinetic equations for the densities of active tips A( with the no-flux boundary condition (∇A = 0, etc.). Here, the duct expansion rate, R elong , and tip reactivation rate, R react , have been replaced by the system expansion rate R exp . As a consistency check, we can see that, in the non-expanding system (R exp = 0), we recover Eq. (S6) by considering delayed tips as belonging to the inactive compartment, i.e., S −→ I.
To proceed, we now consider the influence of the domain expansion on the mean-field equations (S8). Following Crampin et al. [10], where the effect of domain growth was investigated in the context of spatial pattern formation, we considered a general particle density field C(x, t), whose dynamics is governed by a reaction-diffusion equation of the form with R(C) representing the local reaction terms. A domain Ω t , whose topology may depend on time t, introduces a flow field u(x, t), whose effects are captured by a term of the form ∇ · (uC). In this case, the dynamics of the density field C on an expanding domain is governed by the equation [10] ∂C ∂t + ∇ · (uC) = D∇ 2 C + R(C).
Applied to Eqs. (S8), this leads to the mean-field equations The additional terms introduce two new contributions to each field equation: i) advection terms of the form u · ∇C, and ii) a dilution term of the form C(∇ · u), which when expanded leads to Eq. (1) of the main text.
Experimental observations suggest that the mesenchymal tissue undergoes sustained and spatially uniform proliferation, consistent with the observed exponential-like expansion of the developing gland. Here, we assume that growth occurs uniformly across the tissue so that the expansion is affine, preserving locally geometric properties and thus duct morphology. At the same time, we assume that epithelial cells proliferate to preserve the continuity of ducts and branch angles. With this in mind, we may now describe the flow u caused by the expanding domain Ω t as a temporal evolution of the spatial coordinates x(t), so that u = dx dt . Considering a uniform, isotropic, domain expansion, we can relate x(t) at time t, with the coordinates y = x(0) of the original domain Ω 0 by a space-independent function ω(t) [11] x(t) = ω(t)y (S12) where ω(t = 0) = 1, such that For an exponential expansion of the domain, ω(t) ∝ exp(R exp t), with rate constant R exp , we can write the divergence of the flow, or dilution factor, as [11,12] for spatial dimensions d ≥ 1, where R exp d defines the volumetric expansion rate. Note that only for an exponential expansion at constant rate is the ratioω(t)/ω(t) time-independent. Using Eq. (S14), rescaling the rates by the duct creation rate R duct and time by R −1 duct , space by (D/R duct ) 1/2 , defining the dimensionless flow field v = u/(R duct D) 1/2 , and introducing the rescaled density fields a = AR delay /R branch (similarly for S and I), the kinetic equations (1) are written in the dimensionless form as ∂a ∂t +v · ∇a + r exp da = ∇ 2 a + r branch a(1 − (a + i + s)) + r exp s (S15a) ∂s ∂t +v · ∇s + r exp ds = r branch a(a + i + s) − r exp s (S15b) From this result, the effect of the spatial dimension of the system is made clear: as ducts (I) always co-expand with the system, their dilution occurs at a lower rate than that of active (A) and delayed (S) tips. Note that, in dimension d = 1, ducts (I) undergo no dilution at all since, in this case, the ducts and the system expand identically.

Transformation of the mean-field equations to a co-moving frame
Based on this result, we can now represent the reaction-diffusion system (S11) on an isotropically growing domain as a system with time-dependant parameters on a stationary domain [10]. Expressing the active tip density field in terms of the expanding variables a(x, t) = a(ω(t)y, t), and equivalently for the i and s fields, with ω = exp(r exp t), Eqs. (S11) may be written as ∂a ∂t +r exp da = 1 ω 2 (t) ∇ 2 a + r branch a(1 − (a + i + s)) + r exp s (S16a) ∂s ∂t +r exp (d + 1)s = r branch A(a + i + s) (S16b) where the advection term is eliminated in favor of a time-dependent effective diffusion constant that decays as 1/ω 2 (t) = exp(−2r exp t).

Homogeneous solutions
The system of equations (S16) has two homogeneous steady-state solutions: an empty state ρ 0 = (a 0 , s 0 , i 0 ) = (0, 0, 0), and an active state ρ + = (a + , s + , i + ) that exists for d > 1, where . (S17c) This solution is stable when the branching rate is larger than the volumetric expansion rate, i.e., r branch ≥ dr exp (or, equivalently, R branch ≥ dR exp ), and both branching and expansion rates remain non-zero. For the non-expanding case, the only homogeneous steady-state is the empty state ρ 0 , as obtained in the standard BARW model [3]. In Eq. (S16), the empty state ρ 0 becomes unstable when r branch > dr exp , such that any (positive) perturbation away from ρ 0 grows toward the positive homogeneous steady state ρ + . Fronts connecting both homogeneous solutions, ρ 0 and ρ + , may then propagate onto the unstable state, leaving in its wake uniform densities of ducts, and active and inactive tips. The solutions (S17) imply that, for a given dimension d > 1, the ratio of active tips to ducts i + /a + = r exp (d − 1) is controlled by the expansion to duct creation rate ratio r exp , while the ratio of ducts to overall tip density i + /(a + + s + ) = d/(r branch (d − 1)) is inversely proportional to the branching to duct creation rate ratio r branch .

C. Numerical integration of the mean-field equations
The numerical solution of the mean-field equations (S16) was obtained using a finite difference scheme and a 9point stencil for the discrete approximation of the Laplacian (diffusion) operator. For the numerical results shown in Fig. 4l, we considered a lattice spacing dx = 0.1, a time-discretization dt = 0.01, measured in units of (R duct /D) 1/2 and 1/R duct , respectively, a system size of L = 200 points, and dimensionless rates r branch = 1 and r exp = r branch /8. As mentioned in the main text, the mean-field rates do not have a trivial correspondence to the microscopic rates, yet we do expect that the relative ratio r branch /r exp ≈ 8 to be conserved. We considered no-flux boundary conditions at x = 0 and x = L, and periodic boundaries elsewhere, and the system was initialized with a localized perturbation at x = 0 on the active particle field, i.e., a(x, t = 0) = δ(x) and i(x, t = 0) = s(x, t = 0) = 0. In the regimen considered (r branch > r exp d), the active state ρ + given by Eqs. (S17) is stable so that, once initialized as described above, the system converges rapidly toward that state. The front connecting the stable (ρ + ) and unstable (ρ 0 ) states exhibits a peak of density above ρ + that decays on the back of the front with a damped oscillatory profile (Fig. 4l). Once the front has passed, the system reaches a stable stationary state, supported by a constant rate of tip branching and duct creation, which compensate precisely the diluting effect of the expanding domain.

S3. BRANCHING STATISTICS OF INDEPENDENT LOBES OF THE SMG
To explore the dynamics of branching in the salivary gland, we focused our analysis on the first lobe of the SMG, one of its largest components. To explore branching statistics on different lobes, we took advantage of results from an independent investigation of the role of Kras signalling on lineage specification during salivary gland development (see Data and Code Availability). Here, we made use of an inducible genetic labelling system -the Red2Kras confetti mouse model [13] -to study at clonal density the interaction of RFP+ clones active for Kras with wildtype clones carrying a different confetti color. Glands were induced at clonal induction frequency at E13.5, collected and imaged, and their networks manually segmented at E18.5. Clones mutant for Kras showed evidence of a dynamics different from that of wildtype confetti clones. However, at low (clonal) levels of induction, there was no evidence for large-scale changes in the branching characteristics of glands due to the minority of mutant clones. Notably, comparison of the branching statistics of the E18.5 control (confetti) samples with two Kras E18.4 samples from lobes 2 and 4 of the SMG showed that the latter two lobes fell broadly within the range of inter-sample variability of the confetti system ( Supplementary Fig. 12). The fact that different lobes follow the same trends as the confetti data suggests a robust and conserved pattern of branching morphogenesis across the gland, with the IBDRW model describing the large-scale organisation of the salivary gland networks.

S4. BRANCHING DYNAMICS IN OTHER BRANCHING CONTEXTS
Previous studies of branching dynamics of the mouse mammary gland and kidney ductal epithelium suggested that the large-scale organisation of both tissues were predicted well by the standard (non-inflationary) BARW model [3]. In the case of the mammary gland, branching morphogenesis takes place within a preformed mesenchymal tissue known as the fat pad. In this case, expansion is driven by tip-localized progenitors and, as cells leave the tip, they undergo limited proliferation before exiting cell cycle, giving rise to a simple stratified epithelium. In this case, ducts remain largely non-proliferative during the postnatal phase of development. By contrast, in the case of kidney, the embryonic development of the ductal epithelium does not take place in the context of an inactive preformed mesenchyme, but both populations expand cooperatively during development. Moreover, proliferation kinetics measurements show that neither ductal cells nor tip cells in "submerged" endbuds (i.e., endbuds that are positioned away from the periphery of the growing gland) fall out of cycle. It therefore makes sense to question why the statistics of the branched network can be captured by the BARW paradigm.
To address this question, it is helpful to note that the original BARW model is not altogether distinct from the inflationary BARW model. It is instead a limit of the latter model in which the mesenchmye does not contribute to the local expansion of the ducts. In this limit, where the growth rate vanishes, even if submerged tips would have the capacity to continue branching, since there is no local expansion of the surrounding tissue, steric constraints are never alleviated. But what happens when glands "inflate" during the branching process, but the inflation rate is small? Here, it is necessary to question the meaning of "small". Without inflation, the expansion of the ductal network is driven by active tips that become localized at the periphery of the branching network. Tips that become submerged within the network become rapidly inhibited by steric constraints. In this limit, the radial extent of the network expands linearly in time, as predicted by the mean-field dynamics and illustrated in the stochastic simulations of the BARW model. Indeed, such a linear growth characteristic has been demonstrated in the context of the mouse mammary gland development alongside the restriction of proliferation to endbuds at the periphery of the ductal network.
What happens when the inflation rate is non-vanishing? If the expansion of the mesenchymal tissue provides an additional linear-like radial growth dependence, we would expect some contribution from the reactivation of submerged tips as steric constraints become alleviated by the expansion. However, in the long-term limit, such contributions have only a marginal influence on the statistics of branching since such branching events would be infrequent and dwarfed by primary contributions from branching at the periphery of the growing ductal network. By contrast, when the inflation of the mesenchymal tissue leads to a radial expansion that increases exponentially in time, in the longterm limit, the contribution from submerged tips becomes substantial and competitive with peripheral branching, impacting significantly on the branch statistics. To understand why, it is helpful to consider the likelihood that submerged tips can "catch-up" with the surrounding expanding tissue if the expansion is either linear or exponential. In the long term, when the expansion of the tissue is linear, if a tip propagates at a constant speed relative to the expanding tissue, then no matter how small is the ratio of the expansion and propagation speeds, a submerged tip will always eventually catch up with the expanding network (cf. the well-known problem of an "ant on a rubber rope"), which would necessarily cause a temporal arrest in tip proliferation. By contrast, when the expansion is exponential, in the long-term tips would never fully arrest, but their advance would only be delayed by steric influences. This aligns with the results of the long-term stochastic simulations for the salivary gland. However, one can also question what happens at short or intermediate time scales, when the branching networks are limited in size. In this case, the dynamics under a linear or exponential growth characteristic will converge -in other words, an exponential expansion begins as a linear-like expansion, viz.
It takes a characteristic time (t 1/R exp ) before the nonlinearity begins to impact on the growth dynamics. A comparison of the total salivary gland growth (Fig. 3f) is consistent with an exponential-like expansion, which is why we have taken an exponential growth characteristic as a basis to analyse the model. However, for the first lobe of the SMG (Fig. 3f) and also for the kidney (Supplementary Fig. 8a), the expansion does fit well with a linear growth characteristic.
With this in mind, let us consider the results from previous studies of the embryonic development of the mouse kidney. Measurements of the lateral extension of the kidney indicate that, from E11.5 to E19.5 [14], where the dynamics is characterised by branching morphogenesis, the growth characteristic is approximately linear in time ( Supplementary Fig. 8a). Although such behavior is consistent with the dynamics of a tip-driven BARW, one cannot conclude that there is no attendant growth of the underlying mesenchymal tissue. On the contrary, there is evidence that the mesenchyme does continues to grow during this phase, and that the submerged tips and ducts continue to proliferate [15].
Therefore, to investigate quantitatively the effects of a linear-like inflation of the mesenchymal tissue on the statistical organization of the branching network, we performed stochastic simulations of the conventional BARW (i.e., without inflation of the surrounding mesenchyme) and compared the resulting branching topologies to those produced by variants of the IBDRW in which the time-dependence of the inflation is either linear or exponential in time (Figs. 8b-8e). These simulations were performed in a three-dimensional setting, taking a single stalk as the initial condition. In all cases, we considered isotropic branching, setting all parameters of the model to be identical but for the inflation characteristics. The expansion rate for the exponential inflation was set to 0.45 (in units of the rotational diffusion constant) and, for the linear growth, the expansion was proportional to time, with proportionality constant equal to 0.45.
The results of these simulations (Figs. 8b-8e) show that, for the chosen parameters, when the inflation is linear in time, the tree statistics are largely indistinguishable from those of the conventional BARW. This is because the network expansion is dominated by peripheral tip branching, as in the BARW model ( Supplementary Fig. 8f). More generally, if we would increase the rate of linear gland expansion, there would be a transient phase in which submerged tips contribute significantly to the branching organization. However, in general, based solely on the analysis of branching statistics, it is not possible to discriminate between a BARW and IBDRW with a slow linear-like inflationary expansion of the surrounding mesenchyme. Whether the linear-like expansion of the kidney does indeed allow for a slow persistent reactivation of submerged tips, or chemical influences from maturing ducts and/or surrounding mesenchyme inhibit the resumption of branching of submerged tips, remains a subject for future investigation.
These results show that, in slowly growing branching tissues, the branching dynamics may be well-described within the framework of the non-inflationary BARW. However, if the expansion of the surrounding mesenchymal tissue is super-linear or exponential, the reactivation of submerged tips and remodelling of the ductal epithelium may shape the statistical organization of the branching network in a manner that may conform to the IBDRW model.   The location of the root of two submerged branches are highlighted in black, roots of branches with access to the peripheral region of the tissue are highlighted in blue. c Branching tree representation of the ductal network, with the nodes highlighted as in the top panels. Note that submerged tips are still able to proliferate but undergo only a limited degree of branching over the time window, while endbuds belonging to the peripheral branches maintained a highly branching state. All scale bars correspond to 200µm.     Table S1, with a single stalk as initial conditions and a simulations time of 156 hours, to account for the 6.5 days from E14.5 to E18. We note that the tail of the cumulative distributions of subtree sizes (b), which is not captured by any of the models, corresponds to a single large subtree (see large dark blue subtree in Supplementary Fig. 5e) belonging to the sample with the most branching generations (28, compared to 22 and 21 in samples 1 and 3). This sample-to-sample variability could originate from small differences in the organ size and/or developmental time between different mice. Additionally, the long duct in the tail of the distribution of duct lengths (d), correspond to primary ducts of the samples, which are generally long and thick.   [14], exhibit a linear-like growth characteristic, with a fitted linear growth speed of 10µm/h (R 2 = 0.91), with an exponential growth curve producing a much poorer fit to the data (R 2 = 0.84). b-e, Statistics of (b) distribution of duct lengths, (c) endbud fraction, (d) subtree size, and (e) subtree persistence obtained from the standard BARW model and the IBDRW model with either linear or exponential inflation, with parameters as shown in Table S1, considering a single stalk as initial condition, a simulation time equivalent to 300 hours, and uniformly distributed branching angles φ ( Supplementary Fig. 6a), resulting in an isotropic network. The expansion rate for the exponential inflation was set to R exp = 0.05 h −1 while, for linear growth, the expansion was proportional to time, with the same constant of proportionality. f, Realisation of the IBDRW with linear inflation (inset), and a projection of a 2d section along the center of the network, with active tips highlighted in red and inactive tips shown in blue. Note that active tips localized predominantly at the periphery.  Table S1. For each case, only the indicated parameter was varied.  FIG. 12: Branching statistics of lobes 2 and 4 of the mouse SMG. a-d, Statistical distribution of (a) duct lengths, (b) endbud fraction, (c) subtree size, and (d) subtree persistence for the three control samples (all from lobe 1) and n = 2 samples from lobes 2 and 4. Note that, here, data was collected from red2Kras confetti mice induced at (low) clonal density with a small minority of RFP+ labelled cells mutant for Kras (see Supplementary Note).